Initiating
library(oro.dicom)
library(tidyverse)
library(plotly)
library(data.table)
library(htmlwidgets)
library(coreCT)
library(magick)
library(raster)
library(doParallel)
library(oro.dicom)
library(magick)
library(imager)
ref<-1000*data.frame(hmax=(4.64 +0.77),hmin=(4.64 -0.77), mmax=(3.14 + .65), mmin =(3.14- .65 ))
osis<-path.expand("E:/projects/kaggle/osic-pulmonary-fib/osic-pulmonary-fibrosis-progression")
osistrain<-path.expand("E:/projects/kaggle/osic-pulmonary-fib/osic-pulmonary-fibrosis-progression/train")
osistest<-path.expand("E:/projects/kaggle/osic-pulmonary-fib/osic-pulmonary-fibrosis-progression/test")
osisimg<-path.expand("E:/projects/kaggle/osic-pulmonary-fib/img")
setwd(osis)
#invalids<-"ID00011637202177653955184" #https://www.kaggle.com/rashmibanthia/corrupt-files
trd<-read.csv("train.csv") %>% dplyr::rename(Week=Weeks)
#subset(Patient!= invalids)
ted<-read.csv("test.csv") %>% dplyr::rename(Week=Weeks)
invalids<-c("ID00011637202177653955184","ID00052637202186188008618", "ID00128637202219474716089","ID00132637202222178761324")
dd<-list.files(osistrain)
dd <- setdiff(dd,invalids)
*** load dicom img and convert to HU units
x=dd[1]
There were 32 warnings (use warnings() to see them)
read_x<-function(x,thr){
path<-file.path(osistrain,x)
s<-readDICOM(path)
#img<-s$img[[index]]
ct.slope <- unique(extractHeader(s$hdr, "RescaleSlope"))
ct.int <- unique(extractHeader(s$hdr, "RescaleIntercept"))
Ct.n<- extractHeader(s$hdr, "InstanceNumber")
ss <- lapply(s$img, function(x) x*ct.slope + ct.int)
# ff<-function(x){x [ x < -2000 ]<- NA;return(x)}
# trim background
ss<- lapply(ss,function(x){x[x< thr]<-F;return(x)})
#reorder:
fo<-function(i){
index<-which(Ct.n ==i)
return(ss[[index]])
}
so<-lapply(1:length(ss),fo)
return(so)
}
so<-read_x(x,-2000)
**imager
l.im<-lapply(t(so[1:3]),as.cimg)
img<-as.imlist(l.im) ; names(img)<-NULL
plot(img,axes = F)

img<-t(so[[7]])
im<-as.cimg( img )
plot(im)

*** more imager processing
px<-(im < - 600 & im > -860 )
plot(px)

plot(im)
highlight(px)

NA
NA
plot(im)
px.flood(im,78,296,sigma=.21) %>% highlight

pal <- colorRamp(c("black","white"))
fig <- plot_ly(z = ~ so[[7]], colors = pal) %>%
add_surface(
contours = list(
z = list(
show=TRUE,
usecolormap=TRUE,
highlightcolor="#ff0000",
project=list(z=TRUE)
)
)
)
fig <- fig %>% layout(
scene = list(
camera=list(
eye = list(x=1.87, y=0.88, z=-0.64)
)
)
)
fig
NA
** new threshould to focus on lungs
sfail<-read_x(x,-600)
'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2'signed = FALSE' is only valid for integers of sizes 1 and 2
pal <- colorRamp(c("black","white"))
fig <- plot_ly(z = ~ sfail[[7]], colors = pal) %>%
add_surface(
contours = list(
z = list(
show=TRUE,
usecolormap=TRUE,
highlightcolor="#ff0000",
project=list(z=TRUE)
)
)
)
fig <- fig %>% layout(
scene = list(
camera=list(
eye = list(x=1.87, y=0.88, z=-0.64)
)
)
)
fig
It has cut most of the lungs…that is below -800
**** lets start lung segmentation linear model:
d <- as.data.frame(im)
##Subsamble, fit a linear model
m <- sample_n(d,1e4) %>% lm(value ~ x*y,data=.)
##Correct by removing the trend
im.c <- im-predict(m,d)
out <- threshold(im.c)
plot(out)

out<-threshold(im,"auto")
out <- clean(out,3) %>% imager::fill(7)
plot(im)
highlight(out)

NA
NA
im[ out]<- NA
plot(im)

plot_ly(
z = so[[7]],
#colorscale = list(c(0,0.5,1),c("blue", "white", "red")),
colors=pal,
type = "heatmapgl"
)
NA
NA
a<-kmeans(out,2)
colors<- a$centers[a$cluster[1,],]
im<-as.cimg(so[[7]]) %>% imrotate(-90)
detect.edges <- function(im,sigma=1)
{
isoblur(im,sigma) %>% imgradient("xy") %>% enorm %>% imsplit("c") %>% add
}
edges <- detect.edges(im,2) %>% sqrt
plot(edges)

pmap <- 1/(1+edges) #Priority inv. proportional to gradient magnitude
plot(pmap,main="Priority map") #Nice metal plate effect!

NA
NA
seeds <- imfill(dim=dim(pmap)) #Empty image
seeds[170,233,1,1] <- 1 #Background pixel
seeds[36,255,1,1] <- 2 #Foreground pixel
wt1 <- watershed(seeds,pmap)
plot(wt1,main="Watershed segmentation")
the other lung
seeds <- imfill(dim=dim(pmap)) #Empty image
seeds[355,247,1,1] <- 1 #Background pixel
seeds[56,252,1,1] <- 2 #Foreground pixel
wt2 <- watershed(seeds,pmap)
plot(wt2,main="Watershed segmentation")

NA
NA
now, both
wt <- wt1 + wt2
plot(wt)

NA
NA
wt[wt==3]<-T
wt[wt==4]<-F
lungs<-wt * im
plot(lungs)

NA
NA
dx <- imgradient(im,"x")
dy <- imgradient(im,"y")
grad.mag <- sqrt(dx^2+dy^2)
plot(grad.mag,main="Gradient magnitude")
** lets try kmeans
a<-kmeans(im)
plot(im) %>% highlight(a$centers)
LS0tDQp0aXRsZTogIlIgTm90ZWJvb2siDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQpJbml0aWF0aW5nDQoNCmBgYHtyIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9DQpsaWJyYXJ5KG9yby5kaWNvbSkNCmxpYnJhcnkodGlkeXZlcnNlKQ0KbGlicmFyeShwbG90bHkpDQpsaWJyYXJ5KGRhdGEudGFibGUpDQpsaWJyYXJ5KGh0bWx3aWRnZXRzKQ0KbGlicmFyeShjb3JlQ1QpDQpsaWJyYXJ5KG1hZ2ljaykNCmxpYnJhcnkocmFzdGVyKQ0KbGlicmFyeShkb1BhcmFsbGVsKQ0KbGlicmFyeShvcm8uZGljb20pDQpsaWJyYXJ5KG1hZ2ljaykNCmxpYnJhcnkoaW1hZ2VyKQ0KcmVmPC0xMDAwKmRhdGEuZnJhbWUoaG1heD0oNC42NCArMC43NyksaG1pbj0oNC42NCAtMC43NyksIG1tYXg9KDMuMTQgKyAuNjUpLCBtbWluID0oMy4xNC0gLjY1ICkpDQoNCm9zaXM8LXBhdGguZXhwYW5kKCJFOi9wcm9qZWN0cy9rYWdnbGUvb3NpYy1wdWxtb25hcnktZmliL29zaWMtcHVsbW9uYXJ5LWZpYnJvc2lzLXByb2dyZXNzaW9uIikNCm9zaXN0cmFpbjwtcGF0aC5leHBhbmQoIkU6L3Byb2plY3RzL2thZ2dsZS9vc2ljLXB1bG1vbmFyeS1maWIvb3NpYy1wdWxtb25hcnktZmlicm9zaXMtcHJvZ3Jlc3Npb24vdHJhaW4iKQ0Kb3Npc3Rlc3Q8LXBhdGguZXhwYW5kKCJFOi9wcm9qZWN0cy9rYWdnbGUvb3NpYy1wdWxtb25hcnktZmliL29zaWMtcHVsbW9uYXJ5LWZpYnJvc2lzLXByb2dyZXNzaW9uL3Rlc3QiKQ0Kb3Npc2ltZzwtcGF0aC5leHBhbmQoIkU6L3Byb2plY3RzL2thZ2dsZS9vc2ljLXB1bG1vbmFyeS1maWIvaW1nIikNCg0Kc2V0d2Qob3NpcykNCg0KI2ludmFsaWRzPC0iSUQwMDAxMTYzNzIwMjE3NzY1Mzk1NTE4NCIgICAjaHR0cHM6Ly93d3cua2FnZ2xlLmNvbS9yYXNobWliYW50aGlhL2NvcnJ1cHQtZmlsZXMNCg0KdHJkPC1yZWFkLmNzdigidHJhaW4uY3N2IikgICU+JSBkcGx5cjo6cmVuYW1lKFdlZWs9V2Vla3MpDQojc3Vic2V0KFBhdGllbnQhPSBpbnZhbGlkcykNCnRlZDwtcmVhZC5jc3YoInRlc3QuY3N2IikgJT4lIGRwbHlyOjpyZW5hbWUoV2Vlaz1XZWVrcykNCmludmFsaWRzPC1jKCJJRDAwMDExNjM3MjAyMTc3NjUzOTU1MTg0IiwiSUQwMDA1MjYzNzIwMjE4NjE4ODAwODYxOCIsICJJRDAwMTI4NjM3MjAyMjE5NDc0NzE2MDg5IiwiSUQwMDEzMjYzNzIwMjIyMjE3ODc2MTMyNCIpDQoNCmRkPC1saXN0LmZpbGVzKG9zaXN0cmFpbikNCmRkIDwtIHNldGRpZmYoZGQsaW52YWxpZHMpDQpgYGANCg0KDQoNCioqKiBsb2FkIGRpY29tIGltZyBhbmQgY29udmVydCB0byBIVSB1bml0cw0KDQpgYGB7ciBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQ0KeD1kZFsxXQ0KDQpyZWFkX3g8LWZ1bmN0aW9uKHgsdGhyKXsNCiAgcGF0aDwtZmlsZS5wYXRoKG9zaXN0cmFpbix4KQ0KICAgIHM8LXJlYWRESUNPTShwYXRoKQ0KICAgICNpbWc8LXMkaW1nW1tpbmRleF1dIA0KICAgIA0KICAgIGN0LnNsb3BlIDwtIHVuaXF1ZShleHRyYWN0SGVhZGVyKHMkaGRyLCAiUmVzY2FsZVNsb3BlIikpDQogICAgY3QuaW50ICAgPC0gdW5pcXVlKGV4dHJhY3RIZWFkZXIocyRoZHIsICJSZXNjYWxlSW50ZXJjZXB0IikpIA0KICAgIEN0Lm48LSBleHRyYWN0SGVhZGVyKHMkaGRyLCAiSW5zdGFuY2VOdW1iZXIiKQ0KICAgIHNzIDwtIGxhcHBseShzJGltZywgZnVuY3Rpb24oeCkgeCpjdC5zbG9wZSArIGN0LmludCkgDQogICAjIGZmPC1mdW5jdGlvbih4KXt4IFsgeCA8IC0yMDAwIF08LSBOQTtyZXR1cm4oeCl9DQogICAgIyB0cmltIGJhY2tncm91bmQNCiAgICBzczwtICBsYXBwbHkoc3MsZnVuY3Rpb24oeCl7eFt4PCB0aHJdPC1GO3JldHVybih4KX0pDQogICAgI3Jlb3JkZXI6DQogICAgZm88LWZ1bmN0aW9uKGkpew0KICAgICAgaW5kZXg8LXdoaWNoKEN0Lm4gPT1pKQ0KICAgICAgcmV0dXJuKHNzW1tpbmRleF1dKQ0KICAgIH0NCiAgICBzbzwtbGFwcGx5KDE6bGVuZ3RoKHNzKSxmbykNCiAgICByZXR1cm4oc28pDQp9DQpzbzwtcmVhZF94KHgsLTIwMDApDQogICAgDQpgYGANCioqaW1hZ2VyDQoNCmBgYHtyfQ0KbC5pbTwtbGFwcGx5KHQoc29bMTozXSksYXMuY2ltZykNCmltZzwtYXMuaW1saXN0KGwuaW0pIDsgbmFtZXMoaW1nKTwtTlVMTA0KcGxvdChpbWcsYXhlcyA9IEYpDQoNCmBgYA0KDQoNCmBgYHtyfQ0KaW1nPC10KHNvW1s3XV0pDQoNCg0KDQppbTwtYXMuY2ltZyggaW1nICkNCnBsb3QoaW0pDQpgYGANCg0KKioqIG1vcmUgaW1hZ2VyIHByb2Nlc3NpbmcNCg0KDQpgYGB7cn0NCnB4PC0oaW0gPCAtIDYwMCAmIGltID4gLTg2MCApDQoNCnBsb3QocHgpDQpgYGANCg0KDQoNCg0KYGBge3J9DQpwbG90KGltKQ0KaGlnaGxpZ2h0KHB4KQ0KDQoNCmBgYA0KYGBge3J9DQpwbG90KGltKQ0KcHguZmxvb2QoaW0sNzgsMjk2LHNpZ21hPS4yMSkgJT4lIGhpZ2hsaWdodA0KDQpgYGANCg0KDQoNCg0KDQoNCg0KYGBge3J9DQoNCg0KcGFsIDwtIGNvbG9yUmFtcChjKCJibGFjayIsIndoaXRlIikpDQoNCmZpZyA8LSBwbG90X2x5KHogPSB+IHNvW1s3XV0sIGNvbG9ycyA9IHBhbCkgJT4lDQphZGRfc3VyZmFjZSgNCiAgY29udG91cnMgPSBsaXN0KA0KICAgIHogPSBsaXN0KA0KICAgICAgc2hvdz1UUlVFLA0KICAgICAgdXNlY29sb3JtYXA9VFJVRSwNCiAgICAgIGhpZ2hsaWdodGNvbG9yPSIjZmYwMDAwIiwNCiAgICAgIHByb2plY3Q9bGlzdCh6PVRSVUUpDQogICAgKQ0KICApDQopDQpmaWcgPC0gZmlnICU+JSBsYXlvdXQoDQogIHNjZW5lID0gbGlzdCgNCiAgICBjYW1lcmE9bGlzdCgNCiAgICAgIGV5ZSA9IGxpc3QoeD0xLjg3LCB5PTAuODgsIHo9LTAuNjQpDQogICAgKQ0KICApDQopDQoNCmZpZw0KDQpgYGANCg0KKiogbmV3IHRocmVzaG91bGQgdG8gZm9jdXMgb24gbHVuZ3MNCg0KYGBge3J9DQpzZmFpbDwtcmVhZF94KHgsLTYwMCkNCnBhbCA8LSBjb2xvclJhbXAoYygiYmxhY2siLCJ3aGl0ZSIpKQ0KDQpmaWcgPC0gcGxvdF9seSh6ID0gfiBzZmFpbFtbN11dLCBjb2xvcnMgPSBwYWwpICU+JQ0KYWRkX3N1cmZhY2UoDQogIGNvbnRvdXJzID0gbGlzdCgNCiAgICB6ID0gbGlzdCgNCiAgICAgIHNob3c9VFJVRSwNCiAgICAgIHVzZWNvbG9ybWFwPVRSVUUsDQogICAgICBoaWdobGlnaHRjb2xvcj0iI2ZmMDAwMCIsDQogICAgICBwcm9qZWN0PWxpc3Qoej1UUlVFKQ0KICAgICkNCiAgKQ0KKQ0KZmlnIDwtIGZpZyAlPiUgbGF5b3V0KA0KICBzY2VuZSA9IGxpc3QoDQogICAgY2FtZXJhPWxpc3QoDQogICAgICBleWUgPSBsaXN0KHg9MS44NywgeT0wLjg4LCB6PS0wLjY0KQ0KICAgICkNCiAgKQ0KKQ0KDQpmaWcNCmBgYA0KDQpJdCBoYXMgY3V0IG1vc3Qgb2YgdGhlIGx1bmdzLi4udGhhdCBpcyBiZWxvdyAtODAwDQoNCg0KDQogKioqKiBsZXRzIHN0YXJ0IGx1bmcgc2VnbWVudGF0aW9uDQpsaW5lYXIgbW9kZWw6DQoNCmBgYHtyfQ0KZCA8LSBhcy5kYXRhLmZyYW1lKGltKQ0KIyNTdWJzYW1ibGUsIGZpdCBhIGxpbmVhciBtb2RlbA0KbSA8LSBzYW1wbGVfbihkLDFlNCkgJT4lIGxtKHZhbHVlIH4geCp5LGRhdGE9LikgDQojI0NvcnJlY3QgYnkgcmVtb3ZpbmcgdGhlIHRyZW5kDQppbS5jIDwtIGltLXByZWRpY3QobSxkKQ0Kb3V0IDwtIHRocmVzaG9sZChpbS5jKQ0KcGxvdChvdXQpDQoNCmBgYA0KDQoNCg0KYGBge3J9DQogb3V0PC10aHJlc2hvbGQoaW0sImF1dG8iKQ0Kb3V0IDwtIGNsZWFuKG91dCw1KSAlPiUgaW1hZ2VyOjpmaWxsKDcpDQpwbG90KGltKQ0KaGlnaGxpZ2h0KG91dCkNCg0KDQpgYGANCg0KYGBge3J9DQppbW88LWltDQppbW9bIG91dF08LSBOQQ0KcGxvdChpbW8pDQpgYGANCg0KYGBge3J9DQogcGxvdF9seSgNCiAgICAgICAgICAgIHogPSBzb1tbN11dLA0KICAgICAgICAgICAgI2NvbG9yc2NhbGUgPSBsaXN0KGMoMCwwLjUsMSksYygiYmx1ZSIsICJ3aGl0ZSIsICJyZWQiKSksDQogICAgICAgICAgICBjb2xvcnM9cGFsLA0KICAgICAgICAgICAgdHlwZSA9ICJoZWF0bWFwZ2wiDQogICAgICAgICAgICApDQoNCg0KYGBgDQoNCg0KDQpgYGB7cn0NCmE8LWttZWFucyhvdXQsMikNCmNvbG9yczwtIGEkY2VudGVyc1thJGNsdXN0ZXJbMSxdLF0NCg0KDQpgYGANCg0KDQpgYGB7cn0NCmltPC1hcy5jaW1nKHNvW1s3XV0pICU+JSBpbXJvdGF0ZSgtOTApDQoNCg0KZGV0ZWN0LmVkZ2VzIDwtIGZ1bmN0aW9uKGltLHNpZ21hPTEpDQogICAgew0KICAgICAgICBpc29ibHVyKGltLHNpZ21hKSAlPiUgaW1ncmFkaWVudCgieHkiKSAlPiUgZW5vcm0gJT4lIGltc3BsaXQoImMiKSAlPiUgYWRkDQogICAgfQ0KDQplZGdlcyA8LSBkZXRlY3QuZWRnZXMoaW0sMikgJT4lIHNxcnQgDQpwbG90KGVkZ2VzKQ0KcG1hcCA8LSAxLygxK2VkZ2VzKSAjUHJpb3JpdHkgaW52LiBwcm9wb3J0aW9uYWwgdG8gZ3JhZGllbnQgbWFnbml0dWRlDQpwbG90KHBtYXAsbWFpbj0iUHJpb3JpdHkgbWFwIikgI05pY2UgbWV0YWwgcGxhdGUgZWZmZWN0ISANCg0KDQpgYGANCg0KYGBge3J9DQpzZWVkcyA8LSBpbWZpbGwoZGltPWRpbShwbWFwKSkgI0VtcHR5IGltYWdlDQpzZWVkc1sxNzAsMjMzLDEsMV0gPC0gMSAjQmFja2dyb3VuZCBwaXhlbCANCnNlZWRzWzM2LDI1NSwxLDFdIDwtIDIgI0ZvcmVncm91bmQgcGl4ZWwNCg0Kd3QxIDwtIHdhdGVyc2hlZChzZWVkcyxwbWFwKQ0KcGxvdCh3dDEsbWFpbj0iV2F0ZXJzaGVkIHNlZ21lbnRhdGlvbiIpDQoNCg0KDQoNCg0KYGBgDQp0aGUgb3RoZXIgbHVuZw0KDQpgYGB7cn0NCnNlZWRzIDwtIGltZmlsbChkaW09ZGltKHBtYXApKSAjRW1wdHkgaW1hZ2UNCnNlZWRzWzM1NSwyNDcsMSwxXSA8LSAxICNCYWNrZ3JvdW5kIHBpeGVsIA0Kc2VlZHNbNTYsMjUyLDEsMV0gPC0gMiAjRm9yZWdyb3VuZCBwaXhlbA0KDQp3dDIgPC0gd2F0ZXJzaGVkKHNlZWRzLHBtYXApDQpwbG90KHd0MixtYWluPSJXYXRlcnNoZWQgc2VnbWVudGF0aW9uIikNCg0KDQpgYGANCm5vdywgYm90aA0KDQpgYGB7cn0NCg0Kd3QgPC0gd3QxICsgd3QyDQoNCnBsb3Qod3QpDQoNCg0KYGBgDQoNCmBgYHtyfQ0Kd3Rbd3Q9PTNdPC1UDQp3dFt3dD09NF08LUYNCmx1bmdzPC13dCAqIGltDQpwbG90KGx1bmdzKQ0KDQoNCmBgYA0KDQoNCg0KDQoNCg0KDQoNCg0KDQpgYGB7cn0NCg0KZHggPC0gaW1ncmFkaWVudChpbSwieCIpDQpkeSA8LSBpbWdyYWRpZW50KGltLCJ5IikNCmdyYWQubWFnIDwtIHNxcnQoZHheMitkeV4yKQ0KcGxvdChncmFkLm1hZyxtYWluPSJHcmFkaWVudCBtYWduaXR1ZGUiKQ0KDQpgYGANCg0KDQoNCg0KDQoqKiBsZXRzIHRyeSBrbWVhbnMNCg0KDQpgYGB7cn0NCg0KYTwta21lYW5zKGltKQ0KcGxvdChpbSkgJT4lIGhpZ2hsaWdodChhJGNlbnRlcnMpDQoNCg0KYGBgDQpgYGB7cn0NCg0KYGBgDQoNCg==